###############################################################################
# 	File-Name: 	replication
#	Log-file:	na
#	Date:  		05/22/2020
#	Author: 	Jennifer Richmond, Shalu Agrawal, and Johannes Urpelainen
#	Data Used:  Shalu Agrawal; Nidhi Bali; Johannes Urpelainen; Aseem Mahajan; 
#				Daniel Robert Thomas; Sidhartha Vermani; Ryan Kennedy; Smart 
#				Power India; Initiative for Sustainable Energy Policy, 2019, 
#				"Rural Electricity Demand in India (REDI)"
#	Output		na
#	Purpose:   	replication file for SELECT FIGURES in "Drivers of household 
#               appliance usage: evidence from rural India." Energy for 
#               Sustainable Development.
###############################################################################

## Setting working directory for files and loading general packages
setwd("C:/Users/Jen/Dropbox/Rockefeller Analysis/Manuscript_Jen")
library(stargazer)
library(ggplot2)
library(ggsci)
library(ggpubr)
library(reshape)
library(MASS)
library(lfe)
# Load and attach data
mydata<-read.csv("Input/cleandata.csv", header=TRUE)
attach(mydata)

### FIGURE 2: BOX PLOT OF APPLIANCE USAGE HOURS ###

color.object<-element_text(color="gray45")
color.line<-element_line(color="gray45")
# Creating df of vars of interest
energy.usage<-data.frame(v1=new_usage, v2=lighting_usage, v3=cooling_usage,
                         v4=entertain_usage, v5=housekeep_usage)
# Melting it into a long form
eu.df<-melt(energy.usage)
# Plotting boxplot graph
p<-ggplot(eu.df,aes(x=variable, y=value, fill=variable))+ 
    geom_boxplot(alpha=0.25)+
    theme(text = element_text(size=19, color="gray45"))+ 
    theme(axis.text=color.object)+  
    theme(axis.ticks=color.line)+
    scale_x_discrete(labels=c("v1" = "Total", "v2" = "Lighting","v3" = "Cooling", "v4" = "Entertainment", "v5" = "Housekeeping"))+
    ylim(0,100)+labs(y="Hours",x="")+
    guides(fill=FALSE)+
    ggtitle("Totaled hours of household appliance usage by group")
p

### FIGURE 1: Histograms (logged and discrete) of total energy usage ###

# Creating data frame with new_usage var
usage.df<-data.frame(v1=new_usage+1)
#usage.df$v2<-log10(usage.df$v1+1)
usage.df<-melt(v1)
# Creating the hist 
q<-ggplot(usage.df,aes(x=v1, fill=v1))+ 
    geom_histogram(alpha=0.25,fill="#56B4E9", color="#56B4E9")+
    scale_x_log10(breaks=c(0,1,5,25,100,400),labels=c(0,1,5,25,100,400))+
    theme(text = element_text(size=20, color="gray45"))+
    theme(axis.text=color.object)+
    theme(axis.ticks=color.line)+
    labs(y="Households (n=10,249)",x="Total hours/day")+
    ggtitle("Appliance hrs. (ln+1)")
q

usage.df$v1<-usage.df$v1-1
r<-ggplot(usage.df,aes(x=v1, fill=v1))+ 
    geom_histogram(alpha=0.25,fill="#56B4E9", color="#56B4E9")+
    scale_x_continuous(breaks=c(0,100,200,300,400),labels=c(0,100,200,300,400))+
    theme(text = element_text(size=20, color="gray45"))+
    theme(axis.text=color.object)+
    theme(axis.ticks=color.line)+
    labs(y="Households (n=10,249)",x="Total hours/day")+
    ggtitle("Appliance hrs.")
r
figure <- ggarrange(q,r,
                    labels = c("", ""),
                    ncol = 2, nrow = 1)
figure

### FIGURE 6: Coefficients plot ####

# Running OLS models with standardized vars and clustered SEs using lfe package
# 3.1) Lighting 
ols1<-felm(std_light ~ std_grid+std_mg+std_shs+std_backup+std_avail+std_higher_ed+
             std_basic_ed+std_adults+std_children+std_gender+std_hindu+std_caste+
             std_expend+std_dist+std_landowner+std_salaried |q007_state|0|q011_village_code, data=mydata)
summary(ols1)
stargazer(ols1, type="latex")

# 3.2) Cooling
ols2<-felm(std_cool ~ std_grid+std_mg+std_shs+std_backup+std_avail+std_higher_ed+
             std_basic_ed+std_adults+std_children+std_gender+std_hindu+std_caste+
             std_expend+std_dist+std_landowner+std_salaried |q007_state|0|q011_village_code, data=mydata)
summary(ols2)
stargazer(ols2, type="latex")

# 3.3) Entertainment 
ols3<-felm(std_entertain ~ std_grid+std_mg+std_shs+std_backup+std_avail+std_higher_ed+
             std_basic_ed+std_adults+std_children+std_gender+std_hindu+std_caste+
             std_expend+std_dist+std_landowner+std_salaried |q007_state|0|q011_village_code, data=mydata)
summary(ols3)
stargazer(ols3, type="latex")

# 3.4) Housekeeping
ols4<-felm(std_housekeep ~ std_grid+std_mg+std_shs+std_backup+std_avail+std_higher_ed+
             std_basic_ed+std_adults+std_children+std_gender+std_hindu+std_caste+
             std_expend+std_dist+std_landowner+std_salaried |q007_state|0|q011_village_code, data=mydata)
summary(ols4)
stargazer(ols4, type="latex")

# Plotting
library(coefplot)
library(ggplot2)
names<-c("1. Lighting","2. Cooling","3. Entertainment","4. Housekeeping")
color<-c("turquoise4","darkgoldenrod1","coral1","cadetblue3")
multiplot(ols1,ols2,ols3,ols4, single=FALSE, scales="fixed", names=names,
          title="Coefficient plots for each appliance group",
          xlab="Standardized coefficient value", ylab="Variables", decreasing=TRUE,
          legend.position="none", newNames=c(std_grid="Grid",std_mg="Mini-grid",
          std_shs="SHS",std_backup="Backup",std_avail="Electricity availability",std_higher_ed="Higher ed.",
          std_basic_ed="Basic ed.",std_adults="No. adults",std_children="No. children",
          std_gender="Female",std_hindu="Hindu",std_caste="Sched. caste",
          std_expend="Expenditure (ln)",std_dist="Distance to town",
          std_landowner="Landowner",std_salaried="Salaried job"),pointSize=3.7, cex=3, lwdOuter=1)
